En el archivo sat.txt se tienen datos de valores de espectros de pixels en una imagen de satélite utilizados para la predicción de la clase de suelo. Realice un análisis de los datos utilizando diferentes métodos y compárelos mediante Validación Cruzada. Luego aplíquelos a la muestra de test sat.tst y analice los resultados obtenidos.
La base de datos contiene los valores multi-espectrales de píxeles en regiones de 3x3 en una imágen satelital, y la clasificación asociada con el píxel central de cada región. El objetivo es predecir esta clasificación dados los valores multi-espectrales.
Esta base de datos se generó a partir del escáner multi-espectral LANDSAT. Un marco de imágenes Landsat MSS consta de cuatro imágenes digitales de la misma escena en diferentes bandas espectrales. Dos de estos están en la región visible (correspondiente aproximadamente a las regiones verdes y rojas del espectro visible) y dos están en el infrarrojo. Cada píxel es una palabra binaria de 8 bits, con 0 correspondiente al negro y 255 a blanco. La resolución espacial de un píxel es de unos 80 m x 80m. Cada imagen contiene 2340 x 3380 de tales píxeles.
La base de datos es una subárea (pequeña) de una escena, que consta de 82 x 100 píxeles. Cada línea de datos corresponde a una región de píxeles cuadrados 3x3 completamente contenida dentro de la subárea 82x100. Cada línea contiene los valores de píxeles en las cuatro bandas espectrales (convertidas en ASCII) de cada uno de los 9 píxeles en la región 3x3 y un número que indica la etiqueta de clasificación del píxel central.
La clase de cada píxel se codifica como un número, siendo estos:
1: suelo rojo
2: cultivo de algodón
3: suelo gris
4: suelo gris húmdeo
5: suelo con rastrojo de vegetación
6: mezcla de suelos
7: suelo gris muy húmedo
En cada línea de datos, los cuatro valores espectrales para el píxel superior izquierdo se dan primero, seguidos por los cuatro valores espectrales para el píxel medio superior y luego los del píxel superior derecho, y así sucesivamente de modo que los píxeles se leen en secuencia de izquierda a derecha y de arriba a abajo. Por lo tanto, los cuatro valores espectrales para los píxeles centrales están dados por los atributos 17,18,19 y 20. Si lo desea, puede usar solo estos cuatro atributos, ignorando los demás.
Dado que el enunciado menciona que pueden utilizarse solo los datos espectrales correspondientes al píxel central, buscamos extraerlos, junto con su clasificación, a una sola dataframe.
file_trn = "SAT_trn.txt"
data_trn <- read.table(file_trn, sep = "", header=F)
sat_trn_df <- data_trn %>%
dplyr::select(17:20,37) %>%
rename(B1=V17,
B2=V18,
B3=V19,
B4=V20,
C=V37)
file_tst = "SAT_tst.txt"
data_tst <- read.table(file_tst, sep = "", header=F)
sat_tst_df <- data_tst %>%
dplyr::select(17:20,37) %>%
rename(B1=V17,
B2=V18,
B3=V19,
B4=V20,
C=V37)
sat_trn_df$C[sat_trn_df$C == 7] <- 6
sat_tst_df$C[sat_tst_df$C == 7] <- 6
sat_trn_df_numeric <- data.frame(sat_trn_df)
sat_tst_df_numeric <- data.frame(sat_tst_df)
class_names <- c("rojo", "algodón", "gris","gris húmedo", "vegetación", "gris muy húmedo")
class_values <- c(1,2,3,4,5,6)
class_map <- setNames(class_names, class_values)
sat_trn_df$C <- class_map[sat_trn_df$C]
sat_tst_df$C <- class_map[sat_tst_df$C]
| B1 | B2 | B3 | B4 | C |
|---|---|---|---|---|
| 92 | 112 | 118 | 85 | gris |
| 84 | 103 | 104 | 81 | gris |
| 84 | 99 | 104 | 78 | gris |
| 84 | 99 | 104 | 81 | gris |
| 76 | 99 | 104 | 81 | gris |
| 76 | 99 | 108 | 85 | gris |
| 80 | 112 | 118 | 88 | gris |
| 80 | 107 | 113 | 85 | gris |
| 76 | 91 | 104 | 74 | gris húmedo |
| 76 | 95 | 100 | 78 | gris húmedo |
Podemos comprobar que tenemos valores acotados entre 0 y 255 para cuatro bandas espectrales y una categoría de suelo asignada en cada fila. Los contenidos pueden visualizarse mejor utilizando un gráfico de coordenadas paralelas.
Donde cada linea representa una observación para las 4 bandas. La normalización de los valores no es necesaria ya que \(B_1:B_4\) estan en la misma escala. Podemos ver que la distribución de suelo en los datos de entrenamiento y prueba es practicamente la misma.
sat_trn_df <- data.frame(sat_trn_df_numeric)
sat_tst_df <- data.frame(sat_tst_df_numeric)
sat_trn_df$C <- as.factor(sat_trn_df$C)
sat_tst_df$C <- as.factor(sat_tst_df$C)
set.seed(20)
# Entrenamiento del Modelo
n <- nrow(sat_trn_df)
sat_tree <- tree(C ~ B1 + B2 + B3 + B4,
data= sat_trn_df,
control=tree.control(
nobs= nrow(sat_trn_df),
mincut = 1,
minsize = 2,
mindev = 0)
)
# Resultados
sat_tree_summary <- summary(sat_tree)
EEnt <- sat_tree_summary$misclass[1] / n
sat_tree_summary
##
## Classification tree:
## tree(formula = C ~ B1 + B2 + B3 + B4, data = sat_trn_df, control = tree.control(nobs = nrow(sat_trn_df),
## mincut = 1, minsize = 2, mindev = 0))
## Number of terminal nodes: 761
## Residual mean deviance: 0.1723 = 633 / 3674
## Misclassification error rate: 0.04014 = 178 / 4435
sat_tree_CV = cv.tree(sat_tree, FUN = prune.misclass)
size <- sat_tree_CV$size
dev <- sat_tree_CV$dev
sp <-ggplot(data=data.frame(size,dev), aes(x=size, y=dev)) +
geom_line(linetype="solid", color="blue", size=1.2)+
geom_point(color="red", size=3)+
labs(x= "Tamaño del Árbol", y = "Error de CV",title="Poda de Complejidad de Costos")
sp + scale_x_continuous(trans='log10')
# Agrego código para obtener con R el mejor Tamaño
best_size <- which.min(sat_tree_CV$dev)
best_size
## [1] 39
# PQ Harcodea el 13?
sat_tree_prune = prune.misclass(sat_tree, best = 13)
plot(sat_tree_prune)
text(sat_tree_prune, pretty = 0)
CVErr_tree <- sat_tree_CV$dev[39] / n
CVErr_tree
## [1] 0.1815107
set.seed(20)
# Entrenamiento del Modelo con 10-fold CV
sat_rpart <- rpart(C ~ B1 + B2 + B3 + B4, data= sat_trn_df, method = "class", cp=-1,xval=10)
sat_rpart_summary <- printcp(sat_rpart)
##
## Classification tree:
## rpart(formula = C ~ B1 + B2 + B3 + B4, data = sat_trn_df, method = "class",
## cp = -1, xval = 10)
##
## Variables actually used in tree construction:
## [1] B1 B2 B3 B4
##
## Root node error: 3363/4435 = 0.75829
##
## n= 4435
##
## CP nsplit rel error xerror xstd
## 1 2.6167e-01 0 1.00000 1.00000 0.0084779
## 2 2.5335e-01 1 0.73833 0.76896 0.0097636
## 3 1.2905e-01 2 0.48498 0.48498 0.0095487
## 4 6.3634e-02 3 0.35593 0.35742 0.0088020
## 5 1.7544e-02 4 0.29230 0.29379 0.0082400
## 6 1.4570e-02 5 0.27475 0.28189 0.0081181
## 7 6.6905e-03 7 0.24561 0.24234 0.0076694
## 8 5.9471e-03 10 0.22034 0.23253 0.0075467
## 9 5.6497e-03 11 0.21439 0.22629 0.0074660
## 10 4.7577e-03 13 0.20309 0.21885 0.0073673
## 11 4.6090e-03 14 0.19833 0.20934 0.0072364
## 12 3.2709e-03 16 0.18912 0.20636 0.0071945
## 13 1.6354e-03 17 0.18585 0.19774 0.0070698
## 14 1.3381e-03 19 0.18258 0.19625 0.0070479
## 15 1.1894e-03 21 0.17990 0.19536 0.0070346
## 16 1.0903e-03 22 0.17871 0.19477 0.0070258
## 17 8.9206e-04 25 0.17544 0.19715 0.0070610
## 18 6.9382e-04 31 0.17009 0.19715 0.0070610
## 19 5.9471e-04 34 0.16800 0.19715 0.0070610
## 20 4.9559e-04 48 0.15968 0.19685 0.0070567
## 21 4.4603e-04 51 0.15819 0.19387 0.0070124
## 22 3.9647e-04 55 0.15641 0.19328 0.0070035
## 23 2.9735e-04 58 0.15522 0.19328 0.0070035
## 24 1.4868e-04 61 0.15433 0.20012 0.0071047
## 25 8.4958e-05 63 0.15403 0.20101 0.0071176
## 26 7.4338e-05 70 0.15343 0.20071 0.0071133
## 27 0.0000e+00 74 0.15314 0.19952 0.0070960
## 28 -1.0000e+00 197 0.15314 0.19952 0.0070960
plotcp(sat_rpart)
which.min(sat_rpart$cptable[,4])
## 22
## 22
cp<-sat_rpart$cptable[which.min(sat_rpart$cptable[,"xerror"]),"CP"]
cp
## [1] 0.0003964714
CVErr_rpart = sat_rpart$frame[1,'dev'] / n * sat_rpart$cptable[which.min(sat_rpart$cptable[,"xerror"]),"xerror"]
CVErr_rpart
## [1] 0.1465614
poda_rpart<-prune(sat_rpart,cp=cp)
rpart.plot(poda_rpart, cex = 0.5, extra = 0)
set.seed(20)
sat_rf<- randomForest(C~., data= sat_trn_df,
mtry=sqrt(ncol(sat_trn_df)-1), importance=TRUE)
importance(sat_rf)
## 1 2 3 4 5 6
## B1 227.87903 29.96135 362.16076 58.21353 50.21856 48.37936
## B2 140.86987 43.00905 53.51638 37.62347 112.03571 41.93234
## B3 30.04423 10.40657 24.77468 15.11574 28.56435 32.31841
## B4 57.76868 24.80617 34.26451 35.82601 52.16778 74.12139
## MeanDecreaseAccuracy MeanDecreaseGini
## B1 186.14425 1186.7218
## B2 119.11363 914.1809
## B3 48.88156 425.2869
## B4 91.94379 853.1770
varImpPlot(sat_rf)
CVErr_rf <- unname(sat_rf$err.rate[nrow(sat_rf$err.rate),1])
CVErr_rf
## [1] 0.1492672
set.seed(20)
sat_knn <- train(C ~ ., data = sat_trn_df,
method = "knn",
preProcess = c("center","scale"),
trControl = trainControl(method = "cv", number = 10),
tuneGrid = expand.grid(k = 1:100),
metric = "Accuracy")
results <- sat_knn$results
ggplot(results, aes(x = k, y = Accuracy)) +
geom_line(color = "steelblue", size = 1.2) +
xlab("Valor de K") +
ylab("Exactitud") +
ggtitle("Rendimiento para diferentes valores de K") +
theme_bw() +
theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
axis.title.x = element_text(face = "bold", size = 12),
axis.title.y = element_text(face = "bold", size = 12),
axis.text = element_text(size = 10),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank())
best_k <- results$k[which.max(results$Accuracy)]
best_k
## [1] 43
CVErr_KNN <- 1-unname(sat_knn$results[best_k,"Accuracy"])
CVErr_KNN
## [1] 0.1400425
summary(sat_mlr)
## Call:
## nnet::multinom(formula = .outcome ~ ., data = dat, decay = param$decay)
##
## Coefficients:
## (Intercept) B1 B2 B3 B4
## 2 -9.921154 0.5867093 -0.63317495 0.29943995 -0.1122932
## 3 -36.648978 0.7617752 -0.02673514 0.02795903 -0.2534971
## 4 -10.253548 0.5886442 -0.07999703 -0.02004337 -0.2854203
## 5 7.594031 0.4558313 -0.52173889 0.21832973 -0.2099831
## 6 4.599893 0.5678520 -0.18940106 -0.05394639 -0.2887664
##
## Std. Errors:
## (Intercept) B1 B2 B3 B4
## 2 1.6977543 0.05040409 0.04477873 0.05284647 0.04508065
## 3 1.2832527 0.03391253 0.03743092 0.04631104 0.04434272
## 4 0.8182467 0.03437754 0.03467115 0.04409015 0.04040293
## 5 0.9108139 0.03161938 0.03489501 0.04077884 0.03601165
## 6 0.7841282 0.03491517 0.03359844 0.04276092 0.03778148
##
## Residual Deviance: 3755.425
## AIC: 3805.425
#sat_mlr <- nnet::multinom(C ~ ., data = sat_trn_df)
#summary(sat_rl)
# Predict the class labels for the training data
confusion_matrix <- confusionMatrix(predict(sat_mlr, newdata = sat_trn_df), sat_trn_df$C)
CVErr_mlr2 <- 1 - unname(confusion_matrix$overall['Accuracy'])
CVErr_mlr2
## [1] 0.1569335
predicted_mlr <- predict(sat_mlr, newdata = sat_trn_df)
CVErr_mlr <- mean(predicted_mlr != sat_trn_df$C)
CVErr_mlr
## [1] 0.1569335
set.seed(20)
sat_LDA <- stepclass(C ~ B1 + B2 + B3 + B4, data= sat_trn_df, method="lda", criterion = "CR", direction = "both", improvement <- 0.01, fold = 10)
## `stepwise classification', using 10-fold cross-validated correctness rate of method lda'.
## 4435 observations of 4 variables in 6 classes; direction: both
## stop criterion: improvement less than 1%.
## correctness rate: 0.56054; in: "B2"; variables (1): B2
## correctness rate: 0.79776; in: "B1"; variables (2): B2, B1
## correctness rate: 0.82616; in: "B3"; variables (3): B2, B1, B3
##
## hr.elapsed min.elapsed sec.elapsed
## 0.00 0.00 3.13
sat_LDA
## method : lda
## final model : C ~ B1 + B2 + B3
## <environment: 0x000001f8ba1e8210>
##
## correctness rate = 0.8262
CVErr_LDA <- 1-unname(sat_LDA$result.pm["crossval.rate"])
set.seed(20)
sat_QDA <- stepclass(C ~ B1 + B2 + B3 + B4, data= sat_trn_df, method="qda", criterion = "CR", direction = "both", improvement <- 0.01, fold = 10)
## `stepwise classification', using 10-fold cross-validated correctness rate of method qda'.
## 4435 observations of 4 variables in 6 classes; direction: both
## stop criterion: improvement less than 1%.
## correctness rate: 0.58107; in: "B4"; variables (1): B4
## correctness rate: 0.80496; in: "B1"; variables (2): B4, B1
## correctness rate: 0.84849; in: "B2"; variables (3): B4, B1, B2
##
## hr.elapsed min.elapsed sec.elapsed
## 0.00 0.00 3.44
sat_QDA
## method : qda
## final model : C ~ B1 + B2 + B4
## <environment: 0x000001f8bde91e90>
##
## correctness rate = 0.8485
CVErr_QDA <- 1-unname(sat_QDA$result.pm["crossval.rate"])
models <- c("tree", "rpart", "RF", "KNN", "RLM","LDA","QDA")
CVErrs <- c(CVErr_tree,CVErr_rpart,CVErr_rf,CVErr_KNN,CVErr_mlr,CVErr_LDA,CVErr_QDA)
CVErr_df <- data.frame(Label = models, Value = CVErrs)
ggplot(CVErr_df, aes(x = Label, y = Value)) +
geom_bar(stat = "identity", color = "black", fill = "red") + labs(x = "Modelos", y = "Error de Validación Cruzada K=10", title = "Comparación de Precisión de Modelos") +
theme(plot.title = element_text(hjust = 0.5))
#conf_mat <- caret::confusionMatrix(pred_nnet, sat_tst_df$C)
#conf_mat
#conf_matrix(sat_tst_df$C,pred_nnet)
#probs <- predict(model_nnet, newdata = sat_tst_df, type = "probs")
#ROC <- multiclass.roc(response = sat_tst_df$C, predictor = probs)
#auc(ROC)
Queremos comprobar:
sat_trn_df_full <- data_trn
sat_trn_df_full$V37 <- as.factor(sat_trn_df_full$V37)
n <- nrow(sat_trn_df_full)
full_tree <- tree(V37 ~., data= sat_trn_df_full, control=tree.control(nobs= nrow(sat_trn_df_full), mincut = 1, minsize = 2))
full_tree_summary <- summary(full_tree)
EEnt <- full_tree_summary$misclass[1] / n
#full_tree_summary
set.seed(20)
full_tree_CV = cv.tree(full_tree, FUN = prune.misclass)
full_size <- full_tree_CV$size
full_dev <- full_tree_CV$dev
#ggplot(data=data.frame(full_size,full_dev), aes(x=full_size, y=full_dev)) +
# geom_line(linetype="solid", color="blue", size=1.2)+
# geom_point(color="red", size=3)+
# labs(x= "Tamaño del Árbol", y = "Error de CV",title="Poda de Complejidad de Costos")
prune_full = prune.misclass(full_tree, best = 10)
plot(prune_full)
text(prune_full, pretty = 0)
set.seed(20)
full_tree_2 <- rpart(V37 ~., data= sat_trn_df_full, method = "class", cp=-1,xval=10)
full_tree_2_summary <- printcp(full_tree_2)
##
## Classification tree:
## rpart(formula = V37 ~ ., data = sat_trn_df_full, method = "class",
## cp = -1, xval = 10)
##
## Variables actually used in tree construction:
## [1] V1 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V2 V20 V21 V22 V23 V24 V25 V26
## [20] V27 V28 V29 V3 V30 V31 V32 V33 V34 V35 V36 V4 V5 V6 V7 V8 V9
##
## Root node error: 3363/4435 = 0.75829
##
## n= 4435
##
## CP nsplit rel error xerror xstd
## 1 0.26167113 0 1.00000 1.00000 0.0084779
## 2 0.25334523 1 0.73833 0.76896 0.0097636
## 3 0.12905144 2 0.48498 0.48498 0.0095487
## 4 0.06363366 3 0.35593 0.35742 0.0088020
## 5 0.01754386 4 0.29230 0.29497 0.0082518
## 6 0.01486768 5 0.27475 0.28338 0.0081337
## 7 0.01159679 6 0.25989 0.25453 0.0078153
## 8 0.00683913 7 0.24829 0.24591 0.0077128
## 9 0.00669045 9 0.23461 0.23877 0.0076253
## 10 0.00624442 13 0.20458 0.23699 0.0076030
## 11 0.00490633 15 0.19209 0.21588 0.0073269
## 12 0.00356824 17 0.18228 0.21053 0.0072531
## 13 0.00297354 18 0.17871 0.20607 0.0071903
## 14 0.00267618 19 0.17574 0.20517 0.0071776
## 15 0.00237883 21 0.17038 0.20369 0.0071563
## 16 0.00223015 22 0.16800 0.19923 0.0070916
## 17 0.00208147 24 0.16354 0.19923 0.0070916
## 18 0.00178412 26 0.15938 0.19744 0.0070654
## 19 0.00168500 27 0.15760 0.19685 0.0070567
## 20 0.00163544 30 0.15254 0.19685 0.0070567
## 21 0.00148677 32 0.14927 0.19506 0.0070302
## 22 0.00133809 36 0.14332 0.19387 0.0070124
## 23 0.00118941 40 0.13797 0.19447 0.0070213
## 24 0.00104074 50 0.12608 0.19387 0.0070124
## 25 0.00089206 55 0.12073 0.19120 0.0069721
## 26 0.00077312 57 0.11894 0.19060 0.0069631
## 27 0.00074338 62 0.11508 0.19090 0.0069676
## 28 0.00059471 64 0.11359 0.19090 0.0069676
## 29 0.00049559 66 0.11240 0.19179 0.0069811
## 30 0.00029735 69 0.11091 0.19269 0.0069946
## 31 0.00019824 74 0.10943 0.19804 0.0070742
## 32 0.00014868 77 0.10883 0.19893 0.0070873
## 33 0.00000000 83 0.10794 0.19982 0.0071003
## 34 -1.00000000 157 0.10794 0.19982 0.0071003
full_tree_2_NP_pred <- predict(full_tree_2, sat_trn_df_full, type = "class")
# Poda
#full_tree_2$cptable
#plotcp(full_tree_2) #grafico
#which.min(full_tree_2$cptable[,4])
cp<-full_tree_2$cptable[which.min(full_tree_2$cptable[,"xerror"]),"CP"]
poda<-prune(full_tree_2,cp=cp)
rpart.plot(poda, cex = 0.5, extra = 0)
set.seed(20)
rf_full<- randomForest(V37~., data= sat_trn_df_full,
mtry=sqrt(ncol(sat_trn_df_full)-1), importance=TRUE)
#importance(rf_full)
varImpPlot(rf_full)
unname(rf_full$err.rate[nrow(rf_full$err.rate),1])
## [1] 0.08568207